C*********************************************************************
C *** PROGRAM READS NCAR FILES
C *** Input files:  31-36 NEWUPD.1-NEWUPD.6
C ***
C *** Output files: 10    subbox.data
C ***               11    box.data
C ***
C *** This program grids the station data
C*********************************************************************
C****
C**** This program interpolates the given station data or their
C**** ANOMALIES with respect to 1951-1980 to a prescribed grid.
C****
C**** Input files:    units 31,32,...,30+INM (see NCARSURF TODISK)
C****
C****  Record 1: M1,INFOI(2),...,INFOI(8),M1L,TITLEI  header record
C****  Record 2: IDATA(M1-->M1L),LT,LN,ID,HT,NAME,M2,M2L  station 1
C****  Record 3: IDATA(M2-->M2L),LT,LN,ID,HT,NAME,M3,M3L  station 2
C****            etc.
C****  Set ITRIM=0 if the input files are not trimmed (i.e. M1L=
C****  M2L=...= INFOI(4) , INFOI(9) and trailing M2L,... missing)
C****
C**** Output files:     units 10 and 11 (regional means)
C****
C****  10: Record 1: INFO(1),...,INFO(8),TITLEO   header record
C****      Record 2: AVG(1 --> MONM),LTS,LTN,LNW,LNE  grid point 1
C****      Record 3: AVG(1 --> MONM),LTS,LTN,LNW,LNE  grid point 2
C****            etc.
C****  11: Record 1: INFOI(1),...,INFOI(8),TITLER     header record
C****      Record 2: AVGR(1-->M0),WTR(1-->M0),NG,LTS,LTN,LNW,LNE box 1
C****      Record 3: AVGR(1-->M0),WTR(1-->M0),NG,LTS,LTN,LNW,LNE box 2
C****            etc.
C****  AVG(1-->MONM) is a full time series, starting at January
C****  of year IYRBEG and ending at December of year IYREND.
C****  IDATA(1) and AVGR(1) refer to Jan of year IYRBG0 which may
C****  be less than IYRBEG, M0=MONM0 is the max. length of an input
C****  time series, and WTR(M) is the area of the part of the region
C****  that contained valid data for month M (in square meters).
C****
C****  NG is the total number of non-missing data on that record.
C****  LTS,LTN is the latitude of the southern,northern edge of the
C****  (sub)box, LNW,LNE the longitude of the western,eastern edge
C****  in hundredths of degrees (all 4-byte integers).
C****
C****  INFO(1),...,INFO(8)  are 4-byte integers,
C****  TITLEO is an 80-byte character string,
C****  AVG,AVGR,WTR are 4-byte reals.
C****      INFO 1 = 1 (indicates that time series are not trimmed)
C****           2 = KQ (quantity flag, see below)
C****           3 = MAVG (time avg flag: 1 - 4 DJF - SON, 5 ANN,
C****                     6 MONTHLY, 7 SEAS, 8 - 19 JAN - DEC  )
C****           4 = MONM  (length of each time record)
C****           5 = MONM+4 (size of data record length)
C****           6 = IYRBEG (first year of each time record)
C****           7 = flag for missing data
C****           8 = flag for precipitation trace
C****  INFO(I)=INFOI(I) for I=2,3,7,8
C****  In the output file missing data are flagged by
C****  the real number  XBAD = FLOAT( INFO(7) )
C****
C**** The spatial averaging is done as follows:
C**** Stations within RCRIT km of the grid point P contribute
C**** to the mean at P with weight  1.- d/1200,  (d = distance
C**** between station and grid point in  km). To remove the station
C**** bias, station data are shifted before combining them with the
C**** current mean. The shift is such that the means over the time
C**** period they have in common remains unchanged (individually
C**** for each month). If that common period is less than 20(NCRIT)
C**** years, the station is disregarded. To decrease that chance,
C**** stations are combined successively in order of the length of
C**** their time record. A final shift then reverses the mean shift
C**** OR (to get anomalies) causes the 1951-1980 mean to become
C**** zero for each month.
C****
C**** Regional means are computed similarly except that the weight
C**** of a grid box with valid data is set to its area.
C**** Separate programs were written to combine regional data in
C**** the same way, but using the weights saved on unit 11.
C****
C?*** Input parameters (# of input files, time period)
      PARAMETER (INM=6,KM0=12,MONM0=KM0*(2007-1880+1),ITRIM=1,KQM=20)
C?*** Out put parameters (output time period, base period)
      PARAMETER (IYREND=2007, IYBASE=1951,LYBASE=1980)
C?*** Grid dimensions
      PARAMETER (NRM=80,NCM=100)
C?*** Earth radius and limits
      PARAMETER (REARTH=6375.,NCRIT=20,RMAX=999.,RMIN=-999.)
C?*** Work array sizes
      PARAMETER (NSTAM=6000,MSIZE=KM0*300000)
      CHARACTER*80 TITLEI,TITLEO
C**** Input arrays
      DIMENSION IDATA(MSIZE),INFO(8),INFOI(8+ITRIM),ITRL(14+ITRIM)
      REAL RDATA(MSIZE)
C     EQUIVALENCE (IDATA(1),RDATA(1))
C**** Station dependent arrays
      DIMENSION I1S(NSTAM+1),MFS(NSTAM),IORD(NSTAM),ISOFI(NSTAM),
     *   LEN(NSTAM),WTI(NSTAM),ID(NSTAM),ID0(NSTAM),nuseid(nstam)
      REAL*4 CSLONS(NSTAM),CSLATS(NSTAM),SNLONS(NSTAM),SNLATS(NSTAM)
      REAL*4 LONS(NSTAM),LATS(NSTAM) ! lon/lat in degrees
C**** Output-related arrays (dim of output array = MONM < MONM0)
      DIMENSION AVG(MONM0,NCM),DNEW(MONM0),WT(MONM0),WTM(12),
     *  AVGR(MONM0),WTR(MONM0),IORDR(NCM),LENC(NCM)
      COMMON/DIAG/NSTCMB,NSTMNS
C?*** SKIPR and output grid dependent arrays (regions and centers)
      LOGICAL SKIPR/.false./,SKIP(INM,NRM)
      INTEGER INFOR(8)
      REAL*4 PI180,SCALE(KQM)/3*.1,9*0.,.1,7*0./,
     *  CSLONC(NCM,NRM),CSLATC(NCM,NRM),SNLONC(NCM,NRM),SNLATC(NCM,NRM)
      COMMON/GRIDC/PI180,XS(NRM),XN(NRM),XE(NRM),XW(NRM),SKIP,
     *  CSLONC,CSLATC,SNLONC,SNLATC,AREA(NCM,NRM),LATLON(4,NCM+1,NRM)
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      PI180=DATAN(1.D0)/45.D0
      TOMETR=4.*180.*PI180 * REARTH**2/8000.
      RCRIT=1200.
      IF(IARGC().GT.1) THEN
        CALL GETARG(2,TITLEI)
        READ(TITLEI,*) ICRIT
        RCRIT=ICRIT
      END IF
      RBYRC=REARTH/RCRIT
CSRG  RBYRC0=9004.32/RCRIT
      CSCRIT=COS(RCRIT/REARTH)
      NOVRLP=NCRIT
      IYRBEG=1880
      IF(IARGC().GT.0) THEN
        CALL GETARG(1,TITLEI)
        READ(TITLEI,*) IYRBEG
      END IF
C****
C**** Read and use the header record of an input file
C****
      READ (31) INFOI,TITLEI
      KQ=INFOI(2)
      SCL=1.
      IF(KQ.LE.KQM) SCL=SCALE(KQ)
      IF(SCL.EQ.0.) THEN
         WRITE(6,*) ' PROGRAM NOT READY FOR QUANTITY ',KQ
         STOP 'INAPPROPRIATE QUANTITY'
      END IF
C     KM = the number of time frames per year
      KM=1
      IF(INFOI(3).EQ.6) KM=12
      IF(INFOI(3).EQ.7) KM=4
      IF(KM.NE.KM0) STOP 'ERROR: CHANGE KM0'
      ML=INFOI(4)
      NYRSIN=INFOI(4)/KM
      MONM=(IYREND-IYRBEG+1)*KM
      IF(MONM0.NE.INFOI(4).OR.MONM0.LT.MONM) THEN
         WRITE(6,'('' SET MONM0 TO '',I5)') MAX(MONM,INFOI(4))
         STOP 'ERROR: MONM0 NOT OK'
      END IF
      IF(IYRBEG.LT.INFOI(6)) THEN
         WRITE(6,'('' INCREASE IYRBEG TO'',I5)') INFOI(6)
         STOP 'ERROR: IYRBEG TOO SMALL'
      END IF
      MFOUT=1+(IYRBEG-INFOI(6))*KM
      NFB=1+IYBASE-INFOI(6)
      NLB=1+LYBASE-INFOI(6)
      MBAD=INFOI(7)
      LAST=INFOI(7)
      XBAD=MBAD
      TRACE=INFOI(8)
C**** Create and write out the header record of the output file
      DO 10 I=1,8
      INFOR(I)=INFOI(I)
   10 INFO(I)=INFOI(I)
      INFO(1)=1
      INFOR(1)=1
      INFO(4)=MONM
      INFO(5)=INFO(4)+7
      INFOR(5)=2*INFOI(4)+1+4
      INFO(6)=IYRBEG
      TITLEO=TITLEI
      TITLEO(20:39)='ANOM (C)  CR     KM '
      WRITE(TITLEO(33:36),'(I4)') ICRIT
      IF(KQ.EQ.2) TITLEO(25:28)='(MM)'
C?    IF(KQ.EQ.3) TITLEO(31:44)='(MB)          '
C?    IF(KQ.EQ.13) TITLEO(32:50)='(PERCENT)          '
C?    IF(IYBASE.GT.0) TITLEO(61:69)='ANOMALIES'
      WRITE (TITLEO(46:57),'(I4,A)') IYRBEG,'-present'
      WRITE(6,'(1X,A80)') TITLEI,TITLEO
      WRITE(6,'('' INFO BOX   '',8I10)') INFOR
      WRITE(6,'('' INFO SUBBOX'',8I10)') INFO
      WRITE (10) INFO,TITLEO
      IF(.NOT.SKIPR) WRITE(11) INFOR,TITLEO(1:39),TITLEI(40:80)
C**** Find grid-dependent arrays (GRIDC-COMMON)
      CALL GRIDEA (NRM,NCM, INM, RBYRC)
C**** Convert areas to units of square-meters
      DO 20 NR=1,NRM
      DO 20 NC=1,NCM
   20 AREA(NC,NR)=AREA(NC,NR)*TOMETR
C****
C**** Loop over NRM large regions
C****
      DO 300 NR=1,NRM
C**** Collect the station data needed for region NR
      IS=0
      I1S(1)=1
      DO 90 IN=1,INM
      IF(SKIP(IN,NR)) GO TO 90
      REWIND 30+IN
      READ (30+IN) INFOI
      MF=INFOI(1)
      IF(ITRIM.GT.0) ML=INFOI(8+ITRIM)
   50 IF(MF.GE.LAST) GO TO 90
      IF(I1S(IS+1)+ML-MF.GT.MSIZE) STOP 'ERROR: MSIZE TOO SMALL'
      CALL SREAD (30+IN,ML+1-MF,IDATA(I1S(IS+1)),ITRL,14+ITRIM)
      XLAT=.1*ITRL(1)
      XLON=.1*ITRL(2)
      MFS(IS+1)=MF
      LENGTH=ML-MF+1
      MF=ITRL(14)
      IF(ITRIM.GT.0) ML=ITRL(14+ITRIM)
      IF(XLAT.LE.XS(NR).OR.XLAT.GE.XN(NR)) GO TO 50
      XLONR=XLON
      IF(XLON.GT.XE(NR)) XLONR=XLON-360.
      IF(XLON.LT.XW(NR)) XLONR=XLON+360.
      IF(XLONR.LE.XW(NR).OR.XLONR.GE.XE(NR)) GO TO 50
      IS=IS+1
      nuseid(is)=0
      IF(IS.GT.NSTAM) STOP 'ERROR: NSTAM TOO SMALL'
      I1S(IS+1)=I1S(IS)+LENGTH
      LATS(IS)=XLAT
      LONS(IS)=XLON
      CSLATS(IS)=COS(XLAT*PI180)
      SNLATS(IS)=SIN(XLAT*PI180)
      CSLONS(IS)=COS(XLON*PI180)
      SNLONS(IS)=SIN(XLON*PI180)
      ID(IS)=ITRL(3)
      GO TO 50
   90 CONTINUE
      NSTA=IS
      MSIZE0=I1S(NSTA+1)
      IF(NSTA.EQ.0) THEN
         NGOOD=0
         DO 100 M=1,MONM0
         WTR(M)=0.
  100    AVGR(M)=XBAD
         DO 110 N=1,NCM
  110    CALL SWRITE(10,MONM,AVGR,LATLON(1,NC,NR),XBAD)
         GO TO 290
      END IF
C**** Convert data to real numbers if necessary (KQ<100)
      IF(KQ.LT.100) THEN
         DO 115 N=1,I1S(NSTA+1)-1
  115    RDATA(N)=IDATA(N)
      END IF
C**** Order the NSTA stations according to length of time record
      DO 120 IS=1,NSTA
C 120 LEN(IS)=I1S(IS+1)-I1S(IS)
      LEN(IS)=0
      DO 120 M=I1S(IS),I1S(IS+1)-1
C?LIM IF(RDATA(M) >RMAX or <RMIN) RDATA(M)=XBAD
      IF(RDATA(M).EQ.XBAD) GO TO 120
      LEN(IS)=LEN(IS)+1
C?*** Change data if necessary (e.g. trace flag for precip)
C?PRC IF(RDATA(M).EQ.TRACE) RDATA(M)=0.
C?LIM IF(RDATA(M).GT.RMAX or <RMIN) RDATA(M)=RMAX or RMIN
  120 CONTINUE
      CALL SORT (IORD,NSTA,LEN)
C****
C**** Find the time series for all NCM centers in region NR
C****
      DO 200 NC=1,NCM
C**** the subbox edges
      Esou=.01*latlon(1,nc,nr)
      Enor=.01*latlon(2,nc,nr)
      Ewst=.01*latlon(3,nc,nr)
      Eest=.01*latlon(4,nc,nr)
C**** Loop over all stations in memory
      IS0=0
      DO 130 N=1,NSTA
      IS=IORD(N)
C**** Find distance between center and station
      CSDBYR=SNLATS(IS)*SNLATC(NC,NR)+CSLATS(IS)*CSLATC(NC,NR)*
     *  (CSLONS(IS)*CSLONC(NC,NR)+SNLONS(IS)*SNLONC(NC,NR))
CSRG  The next line is replaced by 'IF D/RC>1' below in the orig code
      IF(CSDBYR.LE.CSCRIT) THEN
C****   Keep stations in the current subbox even if D>Rcrit
        dbyrc=1
        if(lats(is).gt.Enor) go to 130
        if(lats(is).lt.Esou) go to 130
        if(abs(lats(is)).lt.89..and.lons(is).gt.Eest) go to 130
        if(abs(lats(is)).lt.89..and.lons(is).lt.Ewst) go to 130
      else
        DBYRC=0.
C**** The arc is replaced by the smaller chord (to avoid using ACOS)
        IF(CSDBYR.LT.1.) DBYRC=RBYRC*SQRT(2.*(1.-CSDBYR))
CSRG    IF(CSDBYR.LT.1.) DBYRC=RBYRC0*SQRT(1.-CSDBYR)
CSRG    IF(DBYRC.GT.1.) GO TO 130
      end if
      IS0=IS0+1
      WTI(IS0)=1.-DBYRC
      ISOFI(IS0)=IS
  130 CONTINUE
C**** Set the weight of the stations outside Rcrit to min(insiders)
      wmin=1.
      do is=1,is0
        if(wti(is).gt.0..and.wti(is).lt.wmin) wmin=wti(is)
      end do
      do is=1,is0
        if(wti(is).le.0.) wti(is)=wmin ! *fac (0<fac<1) ???
      end do
C****
C**** Combine the station data
C****
      NSTCMB=0
      NSTMNS=0
      DO 150 M=1,MONM0
      WT(M)=0.
  150 AVG(M,NC)=XBAD
      IF(IS0.EQ.0) THEN
      LENC(NC)=0
         CALL SWRITE(10,MONM,AVG(1,NC),LATLON(1,NC,NR),XBAD)
         WRITE(6,'('' NO STATIONS FOR CENTER '',2I6)') NR,NC
         GO TO 200
      END IF
C**** Start with the station with the longest time record
      IS=ISOFI(1)
      IOFF=MFS(IS)-I1S(IS)
      NSTMNS=LEN(IS)
      NSTCMB=1
      ID0(1)=ID(IS)
      nuseid(is)=nuseid(is)+1
C**** First update of full data and weight arrays
      WMAX=WTI(1)
      DO 155 K=1,KM
      WTM(K)=WTI(1)
  155 BIAS(K)=0.
      DO 160 M=I1S(IS),I1S(IS+1)-1
      AVG(M+IOFF,NC)=RDATA(M)
      IF(RDATA(M).LT.XBAD) WT(M+IOFF)=WTI(1)
  160 CONTINUE
C**** Add in the remaining stations
      DO 190 I=2,IS0
      IS=ISOFI(I)
      IOFF=MFS(IS)-I1S(IS)
C**** Extend the new data into a full series
      DO 170 M=1,MONM0
  170 DNEW(M)=XBAD
      DO 180 M=I1S(IS),I1S(IS+1)-1
  180 DNEW(M+IOFF)=RDATA(M)
      NF1=1+(MFS(IS)-1)/KM
      NL1=1+(I1S(IS+1)-2+IOFF)/KM
C**** Shift new data, then combine them with current mean
C**** Save the average shifts in the array BIAS
      CALL CMBINE (AVG(1,NC),WT, DNEW,NF1,NL1,WTI(I),WTM, KM,
     *  ID(IS),NSM)
      NSTMNS=NSTMNS+NSM
      IF(NSM.EQ.0) GO TO 190
      NSTCMB=NSTCMB+1
      ID0(NSTCMB)=ID(IS)
      nuseid(is)=nuseid(is)+1
      IF(WMAX.LT.WTI(I)) WMAX=WTI(I)
  190 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      IF(NFB.GT.0) CALL TAVG(AVG(1,NC),KM,NYRSIN, NFB,NLB, NR,NC,0.)
C**** Subtract BIAS, then scale and write the result to disk
      LENC(NC)=0
      M=0
      DO 195 IY=1,MONM0/KM
      DO 195 K=1,KM
      M=M+1
      IF(AVG(M,NC).EQ.XBAD) GO TO 195
      AVG(M,NC)=SCL*(AVG(M,NC)-BIAS(K))
      LENC(NC)=LENC(NC)+1
  195 CONTINUE
      CALL SWRITE(10,MONM,AVG(MFOUT,NC),LATLON(1,NC,NR),RCRIT*(1.-WMAX))
CW    WRITE(6,'('' CENTER '',3I6,'' STATIONS USED'')') NR,NC,IS0
      IF(NSTCMB.EQ.0.OR.NSTCMB.GT.8) GO TO 200
CW    WRITE(6,'(19I7)') (ID(ISOFI(I)),I=1,IS0)
CW    WRITE(6,'(19I7)') (LEN0(ISOFI(I)),I=1,IS0)
      LATC=.5*(LATLON(1,NC,NR)+LATLON(2,NC,NR))
      LONC=.5*(LATLON(3,NC,NR)+LATLON(4,NC,NR))
      WRITE(6,'('' LAT,LON,STN-MNTHS,STNS,IDS'',3I6,I8,2X,8I10)')
     *  LATC,LONC,NSTMNS,NSTCMB,(ID0(I),I=1,NSTCMB)
CW    WRITE(6,'(19I7)') (NINT(1.E6*WTI(I)),I=1,IS0)
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVG(MFOUT-1+M,NC)),M=1,MONM)
  200 CONTINUE
      nuse=0
      do n=1,nstam
      if(nuseid(n).gt.0) then
          write(77,*) 'used station ',id(n),nuseid(n),' times'
          nuse=nuse+1
      end if
      end do
      write(77,*) nuse,' stations used for region', nr
C****
C**** Find the regionally averaged time series
C****
      IF(SKIPR) GO TO 300
      CALL SORT (IORDR,NCM,LENC)
      NC=IORDR(1)
      DO 205 K=1,KM
      WTM(K)=AREA(NC,NR)
  205 BIAS(K)=0.
      DO 210 M=1,MONM0
      WTR(M)=0.
      IF(AVG(M,NC).LT.XBAD) WTR(M)=AREA(NC,NR)
  210 AVGR(M)=AVG(M,NC)
C**** Add in the series of the remaining centers in region NR
      DO 220 N=2,NCM
      NC=IORDR(N)
      IF(LENC(N).EQ.0) GO TO 230
  220 CALL CMBINE (AVGR,WTR,AVG(1,NC),1,MONM0/KM,AREA(NC,NR),WTM,KM,
     *  NC,NSM)
  230 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      IF(NFB.GT.0) CALL TAVG (AVGR,KM,NYRSIN, NFB,NLB, NR,0, 0.)
      NGOOD=0
      M=0
      DO 240 IY=1,MONM0/KM
      DO 240 K=1,KM
      M=M+1
      IF(AVGR(M).EQ.XBAD) GO TO 240
      AVGR(M)=AVGR(M)-BIAS(K)
      NGOOD=NGOOD+1
  240 CONTINUE
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVGR(M)),M=1,MONM0)
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(WTR(M)/TOMETR),M=1,MONM0)
  290 IF(.NOT.SKIPR) WRITE(11) AVGR,WTR,NGOOD,(LATLON(I,NCM+1,NR),I=1,4)
  300 WRITE(6,'('' REGION'',3I9,'' STATIONS USED'')') NR,MSIZE0,NSTA
      STOP
      END

      SUBROUTINE GRIDEA (NRM1,NCM1, INM1, RBYRC)
C****
C**** This output grid dependent routine sets up the latitude and
C**** longitude limits for the Rcrit km hull of the large regions
C**** and computes the relevant grid quantities.
C****
C**** Current order of boxes:   1-4   north , west->east ...
C****                            .      to
C****                          76-80  south , west->east
C**** Order of subboxes:        1-10  south , west->east ...
C****                            .      to
C****                         91-100  north , west->east
C****
      PARAMETER (NRM=80,ICM=10,JCM=10, INM=6,ISHORT=0, NCM=ICM*JCM)
      LOGICAL SKIP(INM,NRM)
      REAL*4 PI180,SNLATJ(JCM),CSLATJ(JCM),LT100S(JCM),LT100N(JCM),
     *  CSLONC(NCM,NRM),CSLATC(NCM,NRM),SNLONC(NCM,NRM),SNLATC(NCM,NRM)
      COMMON/GRIDC/PI180,XS(NRM),XN(NRM),XE(NRM),XW(NRM),SKIP,
     *  CSLONC,CSLATC,SNLONC,SNLATC,AREA(NCM,NRM),LATLON(4,NCM+1,NRM)
C****
C**** Sergej's equal area grid
C****
C**** Grid constants for latitude zones
CCCCC REAL BANDS(9)/90.,64.2,44.4,23.6,0.,-23.6,-44.4,-64.2,-90./
C**** We need only the SINEs of the (northern) band edges
      REAL SNNEDG(9)/1.,.9,.7,.4,0.,-.4,-.7,-.9,-1./
C     Input data sets - 1: 90N-60N  2: 60N-30N ... 6: 60S-90S 7:short
      INTEGER INFRST(8)/1,1,2,2,3,4,5,5/,NUMJ(8)/4,8,12,16,16,12,8,4/
      INTEGER INLAST(8)/2,2,3,4,5,5,6,6/
C**** Check grid dimensions
      IF(NRM1.NE.NRM.OR.ICM*JCM.NE.NCM1.OR.INM1.NE.INM) THEN
         WRITE(6,'(6I9)') NRM1,NRM, ICM*JCM,NCM1, INM1,INM
         STOP 'ERROR: GRID SIZES INCONSISTENT'
      END IF
C**** Don't skip any short files if ISHORT=1
      DO 10 NR=1,NRM
      SKIP(INM,NR)=.FALSE.
      DO 10 IN=1,INM-ISHORT
   10 SKIP(IN,NR)=.TRUE.
C**** Find DDLAT to extend each box R km in north and south direction
      DDLAT=1./(PI180*RBYRC)
C**** Loop over all BOXES (large regions)
      NR=0
      DO 50 J=1,8
C**** Find the sin of the latitudes of the centers in band J
      SNN=SNNEDG(J)
      SNS=SNNEDG(J+1)
      DSLATJ=(SNN-SNS)/JCM
      DO 20 JC=1,JCM
      LT100S(JC)=NINT(100./PI180*ASIN(SNS+(JC-1)*DSLATJ))
      LT100N(JC)=NINT(100./PI180*ASIN(SNS+JC*DSLATJ))
      SNLATJ(JC)=SNS+(JC-.5)*DSLATJ
   20 CSLATJ(JC)=SQRT(1.-SNLATJ(JC)**2)
      DO 50 I=1,NUMJ(J)
      NR=NR+1
C**** Complete definition of SKIP-array
      DO 30 IN=INFRST(J),INLAST(J)
   30 SKIP(IN,NR)=.FALSE.
      DLON=360./NUMJ(J)
      XEAST=-180.+I*DLON
      XWEST=XEAST-DLON
C**** Extend each box by half a box in each direction (>Rcrit km ?)
      XW(NR)=XWEST-.5*DLON
      IF(J.EQ.1.OR.J.EQ.8) XW(NR)=-180.
      XE(NR)=XEAST+.5*DLON
      IF(J.EQ.1.OR.J.EQ.8) XE(NR)=180.
      XN(NR)=ASIN(SNN)/PI180+DDLAT
      XS(NR)=ASIN(SNS)/PI180-DDLAT
C**** Find sine and cosine of the center latitudes and longitudes
      DLONC=DLON/ICM
      DO 40 IC=1,ICM
      LN100W=NINT(100.*(XWEST+(IC-1)*DLONC))
      LN100E=NINT(100.*(XWEST+IC*DLONC))
      RLONI=PI180*(XWEST+(IC-.5)*DLONC)
      SNLONI=SIN(RLONI)
      CSLONI=COS(RLONI)
      DO 40 JC=1,JCM
CSRG  ICJC=JC+JCM*(IC-1)
      ICJC=IC+ICM*(JC-1)
      LATLON(1,ICJC,NR)=LT100S(JC)
      LATLON(2,ICJC,NR)=LT100N(JC)
      LATLON(3,ICJC,NR)=LN100W
      LATLON(4,ICJC,NR)=LN100E
      AREA(ICJC,NR)=1.
      SNLONC(ICJC,NR)=SNLONI
      CSLONC(ICJC,NR)=CSLONI
      SNLATC(ICJC,NR)=SNLATJ(JC)
   40 CSLATC(ICJC,NR)=CSLATJ(JC)
      LATLON(1,NCM+1,NR)=LT100S(1)
      LATLON(2,NCM+1,NR)=LT100N(JCM)
      LATLON(3,NCM+1,NR)=NINT(100.*XWEST)
      LATLON(4,NCM+1,NR)=NINT(100.*(XWEST+DLON))
   50 CONTINUE
      RETURN
      END

      SUBROUTINE CMBINE (AVG,WT, DNEW,NF1,NL1,WT1,WTM, KM, ID,NSM)
C****
C**** Bias of new data is removed by subtracting the difference
C**** over the common domain. Then the new data are averaged in.
C****
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      DIMENSION AVG(KM,*),DNEW(KM,*),WT(KM,*),WTM(*),MISSNG(12)
C**** Loop over months or seasons if appropriate
      NSM=0
      MISSED=KM
      DO 50 K=1,KM
      MISSNG(K)=1
C**** Find means over common domain to compute bias
      SUMN=0
      SUM=0
      NCOM=0
      DO 10 N=NF1,NL1
      IF(AVG(K,N).GE.XBAD.OR.DNEW(K,N).GE.XBAD) GO TO 10
      NCOM=NCOM+1
      SUM=SUM+AVG(K,N)
      SUMN=SUMN+DNEW(K,N)
   10 CONTINUE
      IF(NCOM.LT.NOVRLP) GO TO 50
      BIASK=(SUM-SUMN)/FLOAT(NCOM)
C**** Find mean bias
      WTMNEW=WTM(K)+WT1
      BIAS(K)=(WTM(K)*BIAS(K)+WT1*BIASK)/WTMNEW
      WTM(K)=WTMNEW
C**** Update period of valid data, averages and weights
      DO 20 N=NF1,NL1
      IF(DNEW(K,N).GE.XBAD) GO TO 20
      WTNEW=WT(K,N)+WT1
      AVG(K,N)=(WT(K,N)*AVG(K,N)+WT1*(DNEW(K,N)+BIASK))/WTNEW
      WT(K,N)=WTNEW
      NSM=NSM+1
   20 CONTINUE
      MISSED=MISSED-1
      MISSNG(K)=0
   50 CONTINUE
CLOG? IF(MISSED.GT.0) WRITE(6,90) ID,WT1,MISSNG
   90 FORMAT(' UNUSED DATA - ID/SUBBOX,WT',I8,F5.2,12I2)
      RETURN
      END

      SUBROUTINE TAVG (DATA,KM,NYRS, NFB,NLB, NR,NC, DEFLT)
C****
C**** TAVG computes the time averages (separately for each calendar
C**** month if KM=12) over the base period (year NFB to NLB) and
C**** saves them in BIAS. In case of no data, the average is set to
C**** DEFLT if NR=0 or computed over the whole period if NR>0.
C****
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      DIMENSION DATA(KM,*),LEN(12)
      MISSED=KM
      DO 50 K=1,KM
      BIAS(K)=DEFLT
      SUM=0.
      M=0
      DO 10 N=NFB,NLB
      IF(DATA(K,N).GE.XBAD) GO TO 10
      M=M+1
      SUM=SUM+DATA(K,N)
   10 CONTINUE
      LEN(K)=M
      IF(M.EQ.0) GO TO 50
      BIAS(K)=SUM/FLOAT(M)
      MISSED=MISSED-1
   50 CONTINUE
      IF(NR*MISSED.EQ.0) RETURN
C**** If base period is data free, use bias with respect to whole series
      DO 100 K=1,KM
      IF(LEN(K).GT.0) GO TO 100
      WRITE(6,'(''0NO DATA IN BASE PERIOD - MONTH,NR,NC'',3I9)') K,NR,NC
      SUM=0.
      M=0
      DO 60 N=1,NYRS
      IF(DATA(K,N).GE.XBAD) GO TO 60
      M=M+1
      SUM=SUM+DATA(K,N)
   60 CONTINUE
      IF(M.EQ.0) GO TO 100
      BIAS(K)=SUM/FLOAT(M)
  100 CONTINUE
      RETURN
      END

      SUBROUTINE SORT (INDEX,NDIM,LNGTH)
C**** Sorts INDEX and LNGTH such that LNGTH becomes decreasing
      DIMENSION INDEX(NDIM),LNGTH(NDIM)
      DO 10 N=1,NDIM
   10 INDEX(N)=N
      DO 30 N=1,NDIM-1
C**** Find maximum of LNGTH(N),...LNGTH(NDIM)
      NLMAX=N
      DO 20 NN=N+1,NDIM
      IF(LNGTH(NN).GT.LNGTH(NLMAX)) NLMAX=NN
   20 CONTINUE
C**** Switch positions N and NLMAX
      LMAX=LNGTH(NLMAX)
      LNGTH(NLMAX)=LNGTH(N)
      LNGTH(N)=LMAX
      IMAX=INDEX(NLMAX)
      INDEX(NLMAX)=INDEX(N)
   30 INDEX(N)=IMAX
      RETURN
      END

      SUBROUTINE SREAD (IN,LEN,IDATA,ITRL,NTRL)
C**** Speed read routine for input records
      DIMENSION IDATA(LEN),ITRL(NTRL)
      READ (IN) IDATA,ITRL
      RETURN
      END

      SUBROUTINE SWRITE (IOUT,NDIM,ARRAY,LATLON,DMIN)
      COMMON/DIAG/NSTNS,NSTMNS
      DIMENSION ARRAY(NDIM),LATLON(4)
      WRITE(IOUT) ARRAY,LATLON,NSTNS,NSTMNS,DMIN
      RETURN
      END

      SUBROUTINE LNFIT(DATA,NS,BAD, TREND, IGFRST,IGLAST,NGOOD)
      REAL*8 F1,F2,AX,A1,A2,A22
C**   Determines linear fit using regression analysis
C**   DATA(I) - time series to be fitted
C**   NS - length of the time series
C**   TREND = slope * total number of years
C**   IGFRST,IGLAST = first and last non-missing position
C**   NGOOD = number of good data (trend BAD if NGOOD<2)
      DIMENSION DATA(NS)
      IGFRST=0
      NGOOD=0
      NGSUM=0
      NGSUM2=0
      F1=0.
      F2=0.
      DO 10 I=1,NS
      IF(DATA(I).EQ.BAD) GO TO 10
      IF(IGFRST.EQ.0) IGFRST=I
      IGLAST=I
      NGOOD=NGOOD+1
      NGSUM=NGSUM+I
      NGSUM2=NGSUM2+I**2
      F1=F1+DATA(I)
      F2=F2+DATA(I)*I
   10 CONTINUE
      TREND=BAD
      IF(NGOOD.LT.2) RETURN
      A1=NGOOD
      A2=NGSUM
      A22=NGSUM2
      AX=A1*A22-A2**2
      TREND=NS*(F2*A1-F1*A2)/AX
C     A=(NGSUM2*F1-F2*NGSUM)/AX
      RETURN
      END

